--- title: Peak slice maps (working on it) keywords: fastai sidebar: home_sidebar summary: "Slicing the cube " description: "Slicing the cube " nb_path: "notebooks/60_peakmaps.ipynb" ---
{% raw %}
{% endraw %} {% raw %}
import skimage.exposure as ske
{% endraw %} {% raw %}
from maxrf4u import DataStack, HotmaxAtlas, get_hotslices, get_peakmaps, multi_plot, plot_peakslices, plot_ptrn  

ds = DataStack('RP-T-1898-A-3689.datastack') 

x_keVs = ds.read('maxrf_energies')
y_max = ds.read('maxrf_maxspectrum') 
cube = ds.read('maxrf_cube', compute=False) # don't load into memory yet (too big) 
imvis_reg_highres = ds.read('imvis_reg_highres') 
imvis_reg = ds.read('imvis_reg') 
imvis_extent = ds.read('imvis_extent')

hma = HotmaxAtlas('RP-T-1898-A-3689.datastack')
{% endraw %}

How about Calcium?

{% raw %}
n = 5

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
ymin, ymax = ax.get_ylim()
ax.set_ylim(-5, ymax)
plot_ptrn('Ca', -2, ax);
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack', verbose=True)
slices and super slices:
0: [442 465 488] [442 522]
1: [ 73  96 119] [ 73 119]
2: [472 497 522] [442 522]
3: [1918 1973 2028] [1884 2028]
4: [1884 1905 1926] [1884 2028]
5: [2087 2111 2135] [2087 2135]
6: [270 291 312] [270 312]
7: [945 963 981] [945 981]
8: [1450 1462 1474] [1450 1474]
9: [2194 2215 2236] [2194 2236]
superedge_idxs 16: [73, 119, 270, 312, 442, 522, 945, 981, 1450, 1474, 1884, 2028, 2087, 2135, 2194, 2236]
Computing 16 super slice edges for all baseline corrections...
[########################################] | 100% Completed | 32.5s
Computing slice 0/9
[########################################] | 100% Completed |  3.8s
Computing slice 1/9
[########################################] | 100% Completed |  3.7s
Computing slice 2/9
[########################################] | 100% Completed |  8.2s
Computing slice 3/9
[########################################] | 100% Completed |  5.3s
Computing slice 4/9
[########################################] | 100% Completed |  5.0s
Computing slice 5/9
[########################################] | 100% Completed |  6.1s
Computing slice 6/9
[########################################] | 100% Completed |  6.0s
Computing slice 7/9
[########################################] | 100% Completed |  6.0s
Computing slice 8/9
[########################################] | 100% Completed |  5.5s
Computing slice 9/9
[########################################] | 100% Completed |  4.4s
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

fig, axs = multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel)
{% endraw %}

Let's check if peak #5[6] is an escape peak for Ca_Ka by calculating the energy shift....

{% raw %}
peak_idxs = hma.peak_idxs_list[5]

x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
1.7459311050874269
{% endraw %}

Indeed!!

{% raw %}
title = axs.flatten()[6].get_title() + '(escape)'
axs.flatten()[6].set_title(title);
{% endraw %}

There is much more to be learned from this overview, but before going into this let's look at other elements.

How about Manganese?

Let's create a similar peak map overview for the hotmax spectrum corresponding to Fe_Ka. This is hotmax spectrum #10. Also of interest are the related hotmax spectra for manganese #9 and for Fe_Kb #11.

{% raw %}
n = 9 # Mn

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
[########################################] | 100% Completed | 43.5s
Computing slice 0/10
[########################################] | 100% Completed |  5.6s
Computing slice 1/10
[########################################] | 100% Completed |  5.2s
Computing slice 2/10
[########################################] | 100% Completed | 11.6s
Computing slice 3/10
[########################################] | 100% Completed |  6.4s
Computing slice 4/10
[########################################] | 100% Completed |  5.8s
Computing slice 5/10
[########################################] | 100% Completed |  5.7s
Computing slice 6/10
[########################################] | 100% Completed |  6.1s
Computing slice 7/10
[########################################] | 100% Completed |  6.1s
Computing slice 8/10
[########################################] | 100% Completed |  6.4s
Computing slice 9/10
[########################################] | 100% Completed |  5.8s
Computing slice 10/10
[########################################] | 100% Completed |  5.9s
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

Mm, these histogram equalized peak maps aren't very informative. Most of the slices are just noise. I can not be sure how manganese and iron are correlated. Need to plot that in a different way.

{% raw %}
FeKa_map = peak_maps[3]
MnKa_map = peak_maps[0]
{% endraw %} {% raw %}
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5])

ax.imshow(MnKa_map)
ax.set_title('Mn_Ka')

y, x, z = hot_pixel
ax.scatter([x], [y], c='r', marker='+')
ax1.imshow(FeKa_map, vmax=2)
ax1.set_title('Fe_Ka')
ax1.scatter([x], [y], c='r', marker='+');
{% endraw %}

My first question is, are there any other manganese speckles? And second, do the correlate with iron speckles, as I would expect?

This question needs to wait. Need speckle segmentation functions.

How about Iron?

{% raw %}
n = 10 # Fe_Ka

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
[########################################] | 100% Completed | 29.6s
Computing slice 0/6
[########################################] | 100% Completed |  4.2s
Computing slice 1/6
[########################################] | 100% Completed |  3.2s
Computing slice 2/6
[########################################] | 100% Completed |  5.3s
Computing slice 3/6
[########################################] | 100% Completed |  6.0s
Computing slice 4/6
[########################################] | 100% Completed |  5.6s
Computing slice 5/6
[########################################] | 100% Completed |  6.1s
Computing slice 6/6
[########################################] | 100% Completed |  5.3s
{% endraw %} {% raw %}
ax, labels = hma.plot_spectrum(n) 
plot_peakslices(x_keVs, slices, ax=ax);
ax.set_ylim(-10, 1.15*y_hot.max())
plot_ptrn('Fe', -2, ax);
plot_ptrn('Ca', -1, ax);
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %} {% raw %}
peak_idxs = hma.peak_idxs_list[10]

x_keVs[peak_idxs[0]] - x_keVs[peak_idxs[6]]
1.7459311050874273
{% endraw %}

Again we see an escape peak for Fe_Ka.

How about Chlorine?

{% raw %}
n = 3

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
[########################################] | 100% Completed | 27.0s
Computing slice 0/10
[########################################] | 100% Completed |  3.2s
Computing slice 1/10
[########################################] | 100% Completed |  4.9s
Computing slice 2/10
[########################################] | 100% Completed | 12.3s
Computing slice 3/10
[########################################] | 100% Completed |  6.8s
Computing slice 4/10
[########################################] | 100% Completed |  6.2s
Computing slice 5/10
[########################################] | 100% Completed |  6.4s
Computing slice 6/10
[########################################] | 100% Completed |  6.6s
Computing slice 7/10
[########################################] | 100% Completed |  9.0s
Computing slice 8/10
[########################################] | 100% Completed |  3.0s
Computing slice 9/10
[########################################] | 100% Completed |  6.1s
Computing slice 10/10
[########################################] | 100% Completed |  5.4s
{% endraw %}

We need to make sure we do not confuse Cl_Ka with a Rhodium L peak.

{% raw %}
import moseley as mos 
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax) 
plot_peakslices(x_keVs, slices, ax=ax);

ax1.set_ylim([0, 1.7])

mos.XFluo('Cl', tube_keV=30).plot(ax=ax1, color='b', peak_labels='full')
mos.XFluo('Rh', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
mos.XFluo('K', tube_keV=30).plot(ax=ax1, color='g', peak_labels='simple')
{% endraw %}

This is not the case.

{% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')
{% endraw %} {% raw %}
peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

Seems to be a single grain of salt?

How about Lead and Sulphur?

Run into error for n=16

{% raw %}
n = 18 # Pb 

x_keVs = hma.x_keVs
y_hot = hma.hotmax_spectra[n]
peak_idxs = hma.peak_idxs_list[n] 
hot_pixel =hma.hotmax_pixels[n]

slices = get_hotslices("RP-T-1898-A-3689.datastack", n, clip_level=0.02)
{% endraw %} {% raw %}
peak_maps, keV_maps = get_peakmaps(slices, 'RP-T-1898-A-3689.datastack')
Computing slice 0/7
[########################################] | 100% Completed |  5.8s
Computing slice 1/7
[########################################] | 100% Completed |  3.5s
Computing slice 2/7
[########################################] | 100% Completed |  5.4s
Computing slice 3/7
[########################################] | 100% Completed |  5.9s
Computing slice 4/7
[########################################] | 100% Completed |  6.5s
Computing slice 5/7
[########################################] | 100% Completed | 10.0s
Computing slice 6/7
[########################################] | 100% Completed |  6.1s
Computing slice 7/7
[########################################] | 100% Completed |  6.2s
{% endraw %} {% raw %}
import moseley as mos 
fig, [ax, ax1] = plt.subplots(nrows=2, sharex=True, figsize=[9, 4])
hma.plot_spectrum(n, ax=ax) 
plot_peakslices(x_keVs, slices, ax=ax);

ax1.set_ylim([0, 1.7])

mos.XFluo('Pb', tube_keV=30).plot(ax=ax1, color='k', peak_labels='full')
mos.XFluo('S', tube_keV=30).plot(ax=ax1, color='r', peak_labels='full')
{% endraw %} {% raw %}
titles = [f'[{i}] {x_keVs[peak_idxs[i]]:0.2f} keV' for i in range(len(slices))]
titles.append('VIS')

peak_maps_histeq = [ske.equalize_hist(pm) for pm in peak_maps]

multi_plot(*peak_maps_histeq, imvis_reg, titles=titles, hot_pixel=hot_pixel);
{% endraw %}

This is most interesting. Let's take a look at the keV map for the overlapping band of sulfur and lead. For sulfur this is the S_KL3 band at 2.311 keV. For lead this is the Pb_M5O2 band and 2.372 keV.

The best map for sulfur without lead is actually the potassium map #4

{% raw %}
mid_keV = (2.311 + 2.371) / 2
{% endraw %} {% raw %}
SKa_keV_map = keV_maps[4]
{% endraw %} {% raw %}
fig, [ax, ax1] = plt.subplots(ncols=2, sharex=True, sharey=True, figsize=[9, 5], constrained_layout=True)

ax.imshow(peak_maps[4])
ax.set_title('Overlapping band peak_maps[4]')
pos = ax1.imshow(SKa_keV_map, vmin=2.2, vmax=2.45)
fig.colorbar(pos, ax=ax1, shrink=0.8)
ax1.set_title('Energy of maximum \nS=2.311 keV vs Pb=2.372 keV');
{% endraw %}

Better use best lead free sulfur map in potassium (n=4)?

API

Ok, this is the new Gaussian code base. For now we do not return the fancy peak shape parameters to deconvolve overlapping bands. Soon this will be necessary. For now let's simply take a look at the peak maps...

{% raw %}

gaussian[source]

gaussian(x, x0, sigma)

Normal distribution around `x0` with standard deviation `sigma`.
{% endraw %} {% raw %}

fit_gaussian[source]

fit_gaussian(x, y, peak_idx, rel_height=0.2, baseline=None)

Fit single gaussian distribution at `rel_height`.

Returns: `y_gauss`, `baseline`
{% endraw %} {% raw %}

plot_peakslices[source]

plot_peakslices(x, slices, ax=None, y=None, labels=False)

Plot peak slice regions
{% endraw %} {% raw %}

get_peakmaps[source]

get_peakmaps(slices, datastack_file, norm=False, verbose=False)

Integrate peak `slices`  into peak maps and keV maps.

Returns: peak_maps, keV_maps
{% endraw %} {% raw %}

multi_plot[source]

multi_plot(*images, hot_pixel=None, titles=None, roi_list=None, axis_off=False, sharex=True, sharey=True, vmin=None, vmax=None, cmap='viridis', fontsize='medium', zoom_xyc=None, zoom_half_wh=[100, 100])

Inspect multiple images simultaneously...

Fold along multiple rows if n > 4
{% endraw %} {% raw %}

get_slices[source]

get_slices(x_keVs, y_hot, peak_idxs, rel_height=0.2, clip_level=0.05, baseline=None)

Based on single Gaussian peak fit for each peak slice.

Returns `peak_slices`
{% endraw %} {% raw %}

get_hotslices[source]

get_hotslices(datastack_file, n, rel_height=0.2, clip_level=0.05)

Get slices for specifically for hotmax spectrum `n` in `datastack_file`.

Using noiseline as peakbase.


Returns: `peak_slices`
{% endraw %} {% raw %}

get_superslices[source]

get_superslices(slices, n_channels)

Returns for each slice the superslice edge indices.
{% endraw %} {% raw %}

get_superedges_dict[source]

get_superedges_dict(cube, superslices, verbose=False)

Compute super edges dict (only once).
{% endraw %} {% raw %}

make_baseline[source]

make_baseline(n, slices, superslices, superedges_dict)

Create baseline correction for slice n.

Returns baseline slice
{% endraw %} {% raw %}
{% endraw %}